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NEW  FACTORIZABLE  DISCRETIZATIONS  FOR  THE  EULER  EQUATIONS 

BORIS  DISKIN*  AND  JAMES  L.  THOMASt 


Abstract.  A  multigrid  method  is  defined  as  having  textbook  multigrid  efficiency  (TME)  if  solutions  to 
the  governing  system  of  equations  are  attained  in  a  computational  work  that  is  a  small  (less  than  10)  multiple 
of  the  operation  count  in  one  target-grid  residual  evaluation.  A  way  to  achieve  TME  for  the  Euler  and  Navier- 
Stokes  equations  is  to  apply  the  distributed  relaxation  method  thereby  separating  the  elliptic  and  hyperbolic 
partitions  of  the  equations.  Design  of  a  distributed  relaxation  scheme  can  be  significantly  simplified  if  the 
target  discretization  possesses  two  properties:  (1)  factorizability  and  (2)  consistent  approximations  for  the 
separate  factors.  The  first  property  implies  that  the  discrete  system  determinant  can  be  represented  as  a 
product  of  discrete  factors,  each  of  them  approximating  a  corresponding  factor  of  the  determinant  of  the 
differential  equations.  The  second  property  requires  that  the  discrete  factors  reflect  the  physical  anisotropies, 
be  stable,  and  be  easily  solvable. 

In  this  paper,  discrete  schemes  for  the  nonconservative  Euler  equations  possessing  properties  (1)  and  (2) 
have  been  derived  and  analyzed.  The  accuracy  of  these  scheme  has  been  tested  for  subsonic  flow  regimes 
and  is  comparable  with  the  accuracy  of  standard  schemes.  TME  has  been  demonstrated  in  solving  fully 
subsonic  quasi-one-dimensional  flow  in  a  convergent /divergent  channel. 

Key  words.  Euler  equations,  textbook  multigrid  efficiency,  distributed  relaxation,  factorizable  schemes 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Eull  multigrid  (EMG)  algorithms  [3,  4,  13,  22,  26,  27]  are  the  fastest  solvers  for 
elliptic  problems.  These  algorithms  can  solve  a  general  discretized  elliptic  problem  to  the  discretization 
accuracy  in  a  computational  work  that  is  a  small  (less  than  10)  multiple  of  the  operation  count  in  one 
target-grid  residual  evaluation.  Such  efficiency  is  known  as  textbook  multigrid  efficiency  (TME)  [5,  6]. 
Extending  TME  to  solutions  of  the  Navier-Stokes  equations  is  a  challenging  task  because  the  Navier-Stokes 
equations  form  a  system  of  coupled  nonlinear  equations  that  is  not  fully  elliptic,  even  for  fully  subsonic  flow, 
but  contains  hyperbolic  partitions.  TME  for  the  Navier-Stokes  simulations  can  be  achieved  if  the  different 
factors  contributing  to  the  system  could  be  separated  and  treated  optimally,  e.g.,  by  multigrid  for  elliptic 
factors  and  by  downstream  marching  for  hyperbolic  factors.  One  of  the  ways  to  separate  the  factors  is  the 
distributed  relaxation  method  proposed  in  [3,  4].  The  general  framework  for  achieving  TME  in  large-scale 
computational  fluid  dynamics  (CED)  applications  has  been  discussed  in  [9,  25]. 

The  major  difficulty  in  efficiently  solving  the  Navier-Stokes  equations  is  encountered  with  the  inviscid 
(Euler)  subset;  thus  we  restrict  ourselves  to  the  Euler  equations  here.  The  approach  to  the  solution  of 
the  Euler  equations  motivating  this  paper  is  based  on  an  EMG  algorithm  with  multigrid  cycles  employing 
distributed  relaxation.  It  is  envisioned  that  the  EMG-1  algorithm  (an  EMG  algorithm  with  one  multigrid 
cycle  per  level)  will  provide  solutions  with  algebraic  error  below  the  level  of  the  discretization  error.  Another 
useful  characteristic  of  the  solution  process  is  the  possibility  to  rapidly  converge  residuals  to  the  machine  zero. 
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research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NASl-97046  while 
the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  Virginia  23681 

t Computational  Modeling  and  Simulation  Branch,  Mail  Stop  128,  NASA  Langley  Research  Center,  Hampton,  VA  23681 
(email:  j  .  1 . thomasQlarc  .nasa.gov). 
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The  latter  property  is  not  necessary  for  achieving  TME,  but  it  is  highly  favored  in  practical  applications. 

The  distributed  relaxation  approach  relies  on  a  principal  linearization  of  the  governing  system  of  non¬ 
linear  equations.  The  principal  linearization  is  derived  from  the  full  Newton  linearization  by  removing  some 
unimportant  (subprincipal)  terms.  The  principal  terms  of  a  linear  scalar  equation  are  the  terms  that  make 
major  contributions  to  the  residual  per  a  unit  change  in  the  solution  variable.  The  principal  terms  thus 
generally  depend  on  the  scale,  or  mesh  size,  of  interest.  For  example,  the  discretized  highest  derivative  terms 
are  principal  on  grids  with  small  enough  mesh  size.  For  a  discretized  system  of  differential  equations,  the 
principal  terms  are  those  that  contribute  to  the  principal  terms  of  the  determinant  of  the  matrix  operator. 

Design  of  a  distributed  relaxation  scheme  for  the  Euler  equations  can  be  significantly  simplified  if  the 
target  discretization  possesses  two  properties: 

(1)  The  principal  linearization  of  the  target  discrete  system  is  factorizable  [4,  5,  6,  19,  20],  i.e.,  the 
discrete  system  determinant  can  be  represented  as  a  product  of  discrete  scalar  factors,  each  of  them 
approximating  a  corresponding  factor  of  the  determinant  of  the  differential  equations. 

(2)  The  obtained  scalar  factor  discretizations  reflect  the  physical  anisotropies  and  are  stable  and  easily 
solvable. 

The  main  subject  of  this  paper  is  derivation  of  new  discrete  schemes  for  the  nonconservative  Euler  equations 
possessing  properties  (1)  and  (2).  Corresponding  conservative  discrete  schemes  and  distributed  relaxation 
for  them  have  been  considered  in  [15]. 

Properties  (1)  and  (2)  are  automatically  obtained  with  staggered- grid  discretizations  for  incompressible 
and  slightly  compressible  flow.  Textbook  efficient  multigrid  solvers  employing  factorizable  staggered-grid  dis¬ 
cretizations  of  the  nonconservative  formulations  and  distributed  relaxation  have  already  been  demonstrated 
for  high-Reynolds-number  viscous  incompressible  [11,  24]  and  subsonic  compressible  [23]  flow  regimes. 

Factorizable  schemes  for  the  conservative  Euler  equations  on  collocated  grids  have  been  derived  and 
implemented  in  [17,  18,  19,  20,  21].  The  multigrid  solvers  in  these  references  employed  Collective  Gauss- 
Seidel  rather  than  distributed  relaxation.  The  subsonic-flow  convergence  rates  observed  in  multigrid  V  cycles 
were  quite  fast  (about  0.3  per  cycle),  far  overcoming  the  theoretical  limit  for  nonfactorizable  schemes,  and 
were  only  slightly  grid  dependent.  However,  these  rates  are  still  not  fast  enough  to  guarantee  convergence 
in  an  FMG-1  algorithm.  The  rates  also  deteriorate  somewhat  in  transonic  and  supersonic  computations  and 
for  grids  with  high  aspect  ratios.  These  facts  emphasize  the  need  to  employ  distributed  relaxation. 

This  paper  explores  new  collocated-grid  schemes  for  the  compressible  Euler  equations  that  satisfy  prop¬ 
erties  (1)  and  (2)  in  multiple  dimensions.  A  typical  difficulty  associated  with  this  type  of  scheme  is  a 
poor  measure  of  /i-ellipticity  (stability)  in  the  discrete  approximation  for  the  full-potential  factor  of  the 
system  determinant.  By  definition  (see  [2,  4,  26]),  a  discrete  scalar  (not  necessarily  elliptic)  operator 
L[u]  possesses  a  good  measure  of  /i-ellipticity,  if  the  absolute  value  of  its  symbol  |T(w)j  is  well  sepa¬ 
rated  from  zero  for  all  high-frequency  Fourier  modes.  The  operator  symbol  is  defined  as  the  operator 
response  on  a  discrete  Fourier  mode:  L[e*("''j)]  =  L(w)e*("''j\  where  j  =  {jx,3y,3z)  are  the  grid  indexes 
and  w  =  <  tt  are  normalized  Fourier  frequencies.  For  elliptic  operators, 

high-frequency  Fourier  modes  are  the  modes  satisfying  max(ja;3;j,  jw^j,  \tOz\)  >  7r/2;  for  nonelliptic  operators, 
high-frequency  Fourier  modes  are  those  oscillating  in  the  characteristic  directions. 

Lack  of  /i-ellipticity  often  implies  inefficient  relaxation  (i.e.,  a  poor  smoothing  factor  for  some  high- 
frequency  error  components)  and  slow  convergence  rates  in  multigrid  cycles.  Several  approaches  to  overcome 
the  difficulty  (mainly  in  applications  to  incompressible  flow  equations)  have  been  proposed  (e.g.,  [1,  10]). 
Some  of  the  approaches  are  associated  with  introduction  of  additional  terms  increasing  the  measure  of 
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/i-ellipticity,  others  with  averaging  spurious  oscillations. 

The  factorizable  schemes  for  multidimensional  compressible  flow  equations  proposed  in  this  paper  in¬ 
clude  a  mechanism  to  improve  the  /i-ellipticity  measure  by  obtaining  any  desired  discretizations  for  the 
full-potential  factor.  TME  with  an  FMG  solver  employing  the  distributed  relaxation  method  has  been 
demonstrated  for  two  schemes  approximating  the  Euler  equations  for  a  quasi-one-dimensional  subsonic  flow 
in  a  convergent /divergent  channel. 

The  paper  is  organized  as  follows:  The  Euler  equations  for  inviscid  compressible  flow  problems  are  de¬ 
fined  in  Section  2.  The  idea  of  distributed  relaxation  is  briefly  explained  in  Section  3.  Section  4  presents  the 
derivation  of  the  new  factorizable  schemes  for  the  Euler  equations.  A  model  problem,  the  one-dimensional 
subsonic  Euler  equations,  is  presented  in  Section  5  together  with  a  comparative  analysis  of  linear  discretiza¬ 
tion  schemes.  Description  of  the  multigrid  solver  ingredients  is  in  Section  6  followed  by  results  of  nonlinear 
numerical  tests  with  a  subsonic  flow  in  a  convergent /divergent  channel  reported  in  Section  7.  Section  8 
contains  concluding  remarks. 

2.  Euler  Equations.  The  steady-state  three-dimensional  Euler  system  of  compressible  flow  equations 


can  be  written  as 

udxU  +  vdyU  +  wdzU  +  ^dxP  =  0, 

udxV  +  vdyV  +  wdzV  +  ^dyp  =  0, 

(2.1)  R-(q)  =  <  udxW  +  vdyW  +  wdzW  +  ^dzP  =  0, 

pc^dxU  +  pc^dyV  +  pc^dzW  +  udxP  +  vdyp  +  wdzP  =  0, 

.  ^dxU  + ^dyV  +  ^dzW +  udxe  +  vdye  +  wdze  =  0. 


where  the  primitive  variables,  q  =  {u,v,w,p,e)^ ,  represent  velocity,  pressure,  and  internal  energy,  and  are 
related  to  the  density,  p,  and  the  speed  of  sound,  c,  through  the  following  equations: 

(2.2)  P  =  (7  -  l)pe, 

(2.3)  =  7p/p, 

where  7  is  the  ratio  of  specific  heats. 

In  an  iterative  (quasi-Newton)  procedure,  the  correction  (5q  =  q"+^  —  q",  where  n  is  an  iteration  counter, 
can  be  computed  from  the  equation 

(2.4)  L<5q=-R(q), 
where  L  is  the  principal  linearization  of  the  operator  R(q).  Thus, 

'  Q  0  0  0  ■ 

0  Q  0  ^dy  0 

(2.5)  L  =  0  0  Q  ^dz  0  , 

pc^dx  pc^dy  pc^dz  Q  0 

.  ^dx  ^dy  ^dz  0  Q  _ 

where  Q  =  udx  +  vdy  +  wdz  =  (u  •  V),  and  the  coefficients  u  =  (u,v,w),  p,  and  are  evaluated  from  the 
approximation  q"  and,  for  the  current  iteration,  are  considered  as  constants  unrelated  to  the  target  primitive 
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variables.  The  determinant  of  the  matrix  operator  L  is 
(2.6)  Q^[Q^-c^A], 

where  A  is  the  Laplace  operator.  The  convection  operators,  Q,  and  the  full-potential  operator,  —  c^A, 
represent  hyperbolic  and  elliptic  partitions  of  the  Euler  equations. 

3.  Distributed  Relaxation.  The  distributed  relaxation  method  for  the  Euler  equations  replaces  (5q 
in  (2.4)  by  M(5w,  where 

■  1  0  0  -ia*  0  ■ 

0  1  0  -ia,  0 

(3.1)  M=  0  0  1  -ia^  0  , 

0  0  0  Q  0 

0  0  0  0  1  _ 

so  that  the  resulting  matrix  L  M  becomes  lower  triangular,  as 

'  Q  0  0  0  o' 

0  Q  0  0  0 

(3.2)  LM=0  0  Q  0  0, 

pc^dx  pc^dy  p(?dz  —  c^A  0 

—dx  —dy  —dz  -—A  Q 

and 

(3.3)  LM(5w  =  -R(q). 

The  main  diagonal  of  L  M  is  composed  of  the  factors  of  the  matrix  L  determinant.  The  distributed  relaxation 
approach  yields  fast  convergence  if  the  constituent  scalar  diagonal  operators  in  L  M  are  solved  with  efficient 
methods. 

An  efficient  solver  for  the  convection  factor,  Q,  can  be  based  on  downstream  marching,  with  additional 
special  procedures  for  recirculating  flows  [11,  12,  16,  28].  The  full-potential  factor,  —  c^A,  is  an  operator  of 
variable  type,  and  its  solution  requires  different  procedures  in  subsonic,  transonic,  and  supersonic  regions.  In 
subsonic  regions,  the  full-potential  operator  is  uniformly  elliptic;  therefore  standard  multigrid  methods  yield 
optimal  efficiency.  When  the  Mach  number  approaches  unity,  the  operator  becomes  increasingly  anisotropic 
and,  because  some  smooth  error  components  cannot  be  approximated  adequately  on  coarse  grids,  classical 
multigrid  methods  severely  degrade.  In  supersonic  regions,  the  full-potential  operator  is  uniformly  hyperbolic 
with  the  stream  direction  serving  as  the  time-like  direction.  In  this  region,  an  efficient  solver  can  be  obtained 
with  a  downstream  marching  method.  However,  downstream  marching  becomes  problematic  when  the 
Mach  number  drops  towards  unity,  because  marching  steps  allowed  by  the  stability  condition  are  too  short. 
Thus,  a  special  procedure  is  required  to  provide  an  efficient  solution  for  transonic  regions.  A  possible 
procedure  [7,  8,  14]  is  based  on  piecewise  semicoarsening  and  some  rules  for  adding  dissipation  at  the  coarse 
grid  levels. 

4.  Discrete  Equations.  Having  in  mind  the  distributed  relaxation  procedure  outlined  in  the  previous 
Section  3,  one  would  like  to  design  a  discretization  for  nonlinear  operator  R(q)  of  (2.1)  that  has  the  discrete 
principal  linearization  operator,  Lh,  satisfying  properties  (1)  and  (2)  listed  in  Section  1.  Eor  nonconservative 
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formulations,  the  discretization  of  the  nonlinear  operator  directly  follows  Lh-  Derivation  of  conservative 
discretization  schemes  corresponding  to  a  given  discrete  principal-linearization  operator  has  been  discussed 
in  [15], 

In  this  section,  we  consider  two  factorizable  discretizations  for  the  matrix  operator  L  of  (2.5):  the  basic 
discretization  ,  and  an  improved  discretization,  Lh-  The  basic  collocated-grid  discretization  of 

the  matrix  operator  L  is  defined  as 


Q’^ 

0 

0 

Igh 

p  ^ 

0 

0 

Q'^ 

0 

-d^ 

P  V 

0 

(4.1) 

Th  _ 

basic 

0 

0 

Q'^ 

p  ^ 

0 

pc^dy 

pc^d^ 

0 

i_dh 

7  y 

7  ^ 

0 

Q' 

where  the  discrete  derivatives,  d^,dy,d^,  in  all  off-diagonal  positions  are  the  second-order  accurate  central- 
differencing  approximation.  All  the  diagonal  terms,  except  in  the  fourth  equation,  are  discretized 
with  the  same  second-order  accurate  upwind  (or  upwind-biased)  discretization  scheme.  In  the  subsonic 
regime  (|u|^  =  <  c^),  the  term  is  discretized  with  a  second-order  accurate  downwind  (or 

downwind-biased)  discretization. 

The  determinant  of  the  matrix  operator  is  given  by 

(4.2)  {Q'^y 

where  is  a  wide  (with  mesh  spacing  2h)  discretization  of  the  Laplace  operator.  The  full-potential 
operator  approximation  appearing  in  the  brackets  has  two  major  drawbacks: 

(1)  The  approximation  is  not  /i-elliptic,  i.e.,  it  admits  spurious  oscillatory  solutions  for  the  discrete 
homogeneous  equation. 

(2)  For  near-sonic  regimes  (Mach  number  M  =  |u|/c  ss  1),  the  discrete  operator  stencil  does  not  reflect 
the  physical  anisotropies  of  the  differential  full-potential  operator.  The  discrete  operator  exhibits 
very  strong  coupling  in  the  streamwise  direction,  while  the  differential  operator  has  strong  coupling 
only  in  the  cross-stream  directions. 

An  improved  discrete  full-potential  operator  can  be  obtained  if  the  discretization  is  changed  to 
Then  the  discrete  full-potential  operator  in  (4.2)  becomes 

(4.3) 

If  the  operator  is  second-order  small  (proportional  to  /i^),  the  overall  second-order  discretization  accuracy 
is  not  compromised.  The  choice  of  A^  used  here  is  A'^  =  (q^  -  c^A^/*),  where 

is  a  desired  approximation  for  the  full-potential  factor.  We  do  not  discuss  in  this  paper  the  optimal 
discretization  for  the  multidimensional  subsonic  full-potential  operator.  Note  only  that  it  is  possible  to 
construct  a  discretization  that  satisfies  the  following  properties: 

(1)  For  subsonic  Mach  numbers,  the  discretization  is  /i-elliptic;  in  the  limit  of  the  zero  Mach  number,  it 
is  dominated  by  the  narrow  (with  mesh  spacing  h)  /i-elliptic  Laplace  operator. 

(2)  For  the  Mach  number  approaching  unity,  the  discretization  correctly  reflects  the  physical  anisotropies 
and  tends  to  the  optimal  discretization  for  the  sonic-flow  full-potential  operator  (see  [7,  8,  14]). 

(3)  For  supersonic  Mach  numbers,  the  discretization  becomes  upwind  (or  upwind-biased)  and  can  be 
solved  by  marching. 
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The  operator  is  a  nonlocal  operator  acting  on  and  can  be  introduced  through  a  new  auxiliary 
variable  and  a  new  discrete  equation  =  "D^p^.  Thus,  the  new  vector  of  discrete  unknowns 

becomes  q**  =  (u^,  e^)^.  The  discrete  operator  is  changed  to  Lh,  such  that 

Q'^  0  0  0  0  ' 

0  0  0  0 

0  0  Q'*  0  ia^  0 

0  0  0  Q'*  -V’^  0 

pc^a^  pc^dy  pc^a^  1  0 

c^Qh  ^8^  0  0  Q\ 


The  corresponding  distribution  matrix,  Mh,  is  defined  as 


1  0  0  0  -ia^  0 

0  10  0  -iaf*  0 

p  p 

0  0  10  -iaj  0 

P  Z 

0  0  0  1  v’^  o’ 

0  0  0  0  Q'^  0 

0  0  0  0  0  1 

so  that  the  resulting  matrix  Lh  Mh  becomes  lower  triangular  as 

Q'*  0  0  0  0  o’ 

0  Q'*  0  0  0  0 

0  0  Q'*  0  0  0 

0  0  0  Q'*  0  o' 

pc^a^  pc^dy  pc^d^  1  0 

— aj  — a,^  — a^  o  — 

'y  X  ^  y  ^  Z  'yp 

5.  One-Dimensional  Model  Problem.  The  set  of  the  quasi-one-dimensional  nonconservative  Euler 
equations  is  given  by 

udxU  +  jd^p  =  0, 

(5.1)  pc^dxU  +  udxP  =  --fpu^, 

{-y  -  l)ed:^u  +  ud^e  =  -{'y-l)eu^, 

where  <7(2;)  is  the  area  distribution.  The  principal  linearization  of  the  operator  in  (5.1),  in  the  limit  as  h 
tends  to  zero,  is 

/  udx  0  \ 

(5.2)  pc^dx  udx  0  , 

\  (7  -  l)ea;j  0  ud^  j 

in  which  the  coefficients  u,p,c,  and  e  are  constants  unrelated  to  the  unknown  functions  {u,p,e).  The  third 
equation  is  decoupled  from  the  other  equations.  Thus,  for  the  purpose  of  analysis,  one  can  focus  on  the 
system  of  two  constant-coefficient  equations 

(5.3)  L  q  =  f. 
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where 


(5.4) 


ud^  ^ 

pc^dx  udx  j 


f  =  (/i,  /2)^,  and  q  =  {u,p)'^ .  For  the  subsonic  flow  regimes,  a  natural  set  of  boundary  conditions  for  this 
model  problem  is  u  specified  at  the  inflow  boundary  and  p  specified  at  the  outflow  boundary.  With  this  set 
of  boundary  conditions,  the  differential  problem  (5.3)  is  well  posed. 

The  analysis  presented  in  this  section  compares  the  exact  differential  and  discrete  solutions  for  u  and  p 
obtained  for  the  model  problem  (5.3).  Let  the  exact  solution  of  (5.3)  defined  on  the  interval  x  £  [0, 1]  be 


(5.5)  qexact(a;) 

where  a  is  an  arbitrary  frequency.  Then 


‘^exact 

Pexact 


Cu 

c„ 


iax 


(5.6) 


\  h  J  \  +  uCp  J 


The  system  (5.3)  is  subject  to  boundary  conditions 


(5.7) 

The  distribution  matrix  M, 

(5.8) 

results  in 

(5.9) 


i{0)  =  C^,  p{l)=Cp 


(  1  --dx 
M  =  P 

\  0  udx 


LM  = 


udx  0  \ 

pc^dx  T  J 


The  main  diagonal  of  matrix  LM  is  composed  of  the  convection  operator  udx  and  the  full-potential  operator 
T  =  {u?  —  c^)dxx-  The  one-dimensional  problem  is  very  specific  for  at  least  two  reasons: 

(1)  The  full-potential  factor  vanishes  at  the  sonic  speed  {u  =  c). 

(2)  The  characteristics  perfectly  align  with  the  grid. 

Both  these  features  disappear  in  multiple  dimensions. 

The  corresponding  discrete  problem  is  defined  on  a  uniform  grid  with  mesh  size  h  as 


(5.10) 


-‘h  q 


=  r 


where  Lh  is  a  discretization  of  L,  and  q^  =  {uj,Pj)'^  and  fj*  =  {fi{jh),  are  discrete  representations 

of  the  solutions  and  source  functions,  respectively,  and  j  =  0, 1, 2, . . . ,  77,  77  =  l//i.  The  general  solution  to 
(5.10)  can  be  sought  as  a  combination  of  a  particular  solution  and  the  general  solution  to  the  corresponding 
homogeneous  problem 


(5.11)  Lh  q**  =  0. 

A  particular  solution  can  be  found  in  the  form 


(5.12) 


^Dar 


11 

^par 

Ppai 
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(5.13) 


where  to  =  ah  is  a  normalized  frequency,  and  Lh(A)  is  a  generalized  matrix  symbol  of  the  discrete  operator  Lj,. 
The  entries  of  Lh(A)  are  generalized  symbols  of  the  discrete  scalar  operators  composing  Lh  that  are  defined 
as  responses  of  these  operators  on  the  exponent  function  A-^  .  For  example,  the  generalized  symbol,  of 

the  central  second-order  accurate  difference  approximation  to  the  first  derivative,  3'^,  is  =  3'^(A)A'^, 

n\)  =  i(A-i). 

The  general  solution,  of  the  homogeneous  system  of  equations  (5.11)  is  a  combination  of  linearly 

independent  characteristic  solutions  where  Xk  are  the  roots  of  the  characteristic  polynomial 


(5.14) 


det  Lh(A)  =  0. 


(5-15)  Qhom  =  ( 

\  Phom  J  j  k  k 

The  general  solution  of  the  discrete  problem  (5.10)  is 
(5-16)  q**  =  qhom  +  Qpar  =  Y.  + 

k 

Parameters  Ck  are  chosen  to  satisfy  a  set  of  discrete  boundary  conditions.  The  discretization  error  function 
is  defined  as 

(5-17)  q'^-gexact^ 

where  q^xact  i®  ^  restriction  of  qexact(a;)  to  the  grid  with  the  mesh  size  h. 

Below,  the  discretization  errors  for  four  discrete  factorizable  schemes  approximating  (5.3)  are  compared: 
Scheme  #  1.  The  “basic”  scheme  of  the  (4.1)  type. 

Scheme  #  2.  The  standard  upwind  discrete  scheme. 

Scheme  #  3.  The  discrete  scheme  of  the  (4.4)  type  with  the  discretization  for  the  full-potential 
factor  given  as 

(5.18)  J''*  =  (u^  -  c^)  3^3'^, 

where  the  discrete  operators,  and  3*^,  are  second-order  accurate  upwind  and  downwind  difference 
approximations  to  the  first  derivative,  respectively. 

(4)  Scheme  #  4-  The  discrete  scheme  of  the  (4.4)  type  with  the  discretization  for  the  full-potential 
factor  given  as 

(5.19)  =  (u‘^  -  c‘^) 

where  the  discrete  operator  is  a  three-point  central  approximation  to  the  second  derivative. 

All  the  schemes,  except  the  standard  upwind  scheme  (2),  are  factorizable  in  multiple  dimensions.  The 
discrete  boundary  conditions  for  all  the  four  schemes  are  overspecified,  i.e.,  the  discrete-solution  values  at 
the  boundary  and,  wherever  necessary,  outside  of  the  target  computational  domain  are  specified  from  the 
known  exact  solution  (5.5). 

8 


5.1.  Scheme  #  1.  The  one-dimensional  version  of  the  “basic”  collocated-grid  discretization  for  matrix 
operator  L  of  (5.3)  is  defined  as 

(5.20)  =  f**, 


(5.21) 


gu  Igc 
P 

pc^d^  ud^ 


Recall  that  the  discrete  derivatives,  d^,  and  d'^,  are  second-order  accurate  upwind,  central,  and  downwind 
difference  approximations,  respectively. 

The  generalized  symbol  for  the  operator  is  defined  as 


(5.22) 


where 


ud^{X)  la" (A)  \ 
pc^a"(A)  ud'^{X)  j 


(5.23) 


a^A)  =  i(|-2i  +  i^), 
=  H-f  +  2A-iA2), 
^'(A)  =  H^A-ix). 


A  particular  solution  to  (5.20)  is 

(5.24)  qW  =  (  “11, '  )  e- 


where 

..(1)  _  59'^(e*")/i-i9Te'")/2 

(5.25) 

-(1)  _  -pc^ 9*^(6*"  )/i -1-59"  (e*")/2 


A  set  of  linearly  independent  characteristic  solutions  zi^(j)  =  is  given  by 


(5.26) 

which  corresponds  to  Ai,2  =  1; 

(5.27) 


vi  = 


and  V2  = 


^3  I  piid“(\s] 

\  d-(\s) 


which  corresponds  to  A3  = 


(5.28) 


_  5u^  —  c^—4uVu‘^  —  c‘^  . 


Zu^-\-c^ 


and 


V4  = 


pud^{\4) 

d<^{U) 


which  corresponds  to  A4  =  — — ^^'^2 —  • 

The  characteristic  solutions  z^,  are  normalized  to  satisfy  max|zi;(j)|  =  0(1)  as  h  tends  to  zero.  The 

i 

characteristic  solutions  zi  and  Z2  correspond  to  solutions  of  the  target  differential  problem;  the  characteristic 
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solutions  Z3  and  Z4  are  numerical  artifacts.  Note  that  in  subsonic  regimes,  IA3I  =  IA4I  =  1;  this  implies 
existence  of  global  discrete  solutions  that  do  not  approximate  solutions  of  the  differential  problem.  These 
spurious  solutions  are  a  source  of  instability  of  the  discrete  approximation  (5.20).  Details  are  given  in 
Appendix  A.  For  stable  approximations,  characteristic  solutions  unrelated  to  the  differential  solutions  should 
be  local,  i.e.,  they  should  correspond  to  |Ai;|  1.  The  coefficients  cu  are  computed  from  a  4x4  linear  system 

that  arises  after  substituting  the  general  solution  into  the  boundary  condition  equations 

(i)  Mq  =Cu, 

(ii)  +  =/iW, 

(hi)  -  5<_2)  +  |(-fPN-i =/2(l-/i), 

(iv)  p%  =Cpe*°^, 

where  values  and  p%^i  are  specified  from  the  exact  solution  (5.5). 

5.2.  Scheme  #  2.  A  one-dimensional  version  of  the  standard  upwind  scheme  for  matrix  operator  L 
of  (5.3)  is  defined  as 

(5.30)  =  f**. 


_i_  u  —  cpid  u-\-c  piu  u  —  c  fid 

15  311  +—0  -^0 

'  '  h  pc(u  +  c)  gu  _  pc(u  —  c)  Qd  u+C  QU  _j_  U  —  C  Qd 

The  following  four  boundary  conditions  (two  from  the  left  and  two  from  the  right)  are  used  by  the 
interior  discretizations: 


(i)  ho  + 

(ii)  h-i  +  hcP-i 


—  ir<  4 _ i_r< 

—  ^  2pc^V 

=  h^u  +  6 


(hi)  iu(V+i  -  h^P%+i  =  -  ^Cpj 

(iv)  -  ^^p%  =  -  ^Cp^ 


A  particular  solution  can  be  found  in  the  form 


q(2)  = 

^par 


where 


.(2)  _  + 


p'  '  = 


_  -pcm£d'‘(e‘-)-^d^(e‘-))h  +  m£d'‘(e‘-)  +  ^d^(e‘-))f2 


A  set  of  linearly  independent  characteristic  solutions  zi^(j)  =  v/.Xl  is  given  by 


and  V2  = 


which  corresponds  to  Ai  2  =  1; 
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which  corresponds  to  A3  =  |;  and 
(5.37) 


V4  =  (A4) 


-N 


1 

-pc 


which  corresponds  to  A4  =  3.  The  characteristic  solutions  zi  and  Z2  approximate  solutions  of  the  differential 
equations;  the  characteristic  solutions  Z3  and  Z4  are  local.  The  discrete  scheme  (5.30)  is  stable.  The 
coefficients  ct  are  found  by  substituting  the  general  solution  into  (5.32). 

5.3.  Scheme  #  3.  A  factorizable  scheme  corresponding  to  (4.4)  is  defined  as 


(5.38) 


(5.39) 


(5.40) 


L?>  = 


(  ud^ 

0 


0 

ud^ 

1 


1(3)  = 


=  r 


ud'^  J 


and  the  desired  discrete  full-potential  operator  is  given  by 
(5.41)  =  {u^  -  c^)d^d^. 

The  overspecified  boundary  conditions,  where  values  of  and  p(v_|_i  are 

specified  from  the  exact  differential  solution  (5.5),  and  =  tpo  =  =  0,  are  equivalent  to  the 

following  six  boundary  conditions: 


(5.42) 


In  evaluation  of  (v),  the  value  of  is  computed  from  the  equation  —  'DVjv-i  =  0. 

A  particular  solution  of  the  nonhomogeneous  problem  can  be  found  as 


(i) 

Uq 

=  c^, 

(ii) 

ud^u'i  + 

=  fi(h), 

(iii) 

ud'^ipi  —  V^pi 

=  0, 

(iv) 

=  0, 

(v) 

pc‘^d^u%_^  +  +  ud^P%_i 

=  /2(1  - 

(vi) 

Pn 

=  Cpe-. 

(5.43) 


where 


(5.44) 


/  u(3) 

q<2  =  I  c'“‘, 

\  0' 


_  9‘(e*")  ^-pc^9‘(e*")/i-|-4j9"(e*")/2 


jAZ)  —  ( -pc^9Te*")/i-t449"(e*")/2  ^ 


P  - 
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The  generalized  symbols  d^{X),d^{X),  and  3'^(A)  are  defined  in  (5.23)  and 

=  (u2 -c2)a“(A)a'^(A), 

V^ix)  =  c2  ((a^(A))2  -  a“(A)a'^(A))  =  -  4i  +  e  -  4a  +  a^)  . 

The  six  linearly  independent  characteristic  solutions,  zi^(j)  =  v/^(j)Xl,  are  given  by 


(5.45) 


(5.46) 

Vl  = 

,V2  = 

o  o 

f 

,  and  V3  = 

V 

jh 


\ 


pc‘ 


(#-0 


\  -jhfm  j 


which  corresponds  to  Ai,2,3  =  1, 


/ 

h 

( 

h 

\ 

V4  = 

—p(?hd^{X4) 

j 

= 

4  2 

V 

0 

) 

{ 

0 

y 

(5.47) 


and 


(  jh 

] 

( 

jh 

\ 

V5  = 

pc^h  _  d^{X,)  -  jdHX,)] 

1 

= 

P(p  ( 

( 

1 

-jhpu 

y 

which  corresponds  to  A4,5  =  where  d^(X)  =  ^  (2-^  —  ^)  and  3'=(A)  =  \{W  +  ^A) ; 
and 


/ 


(5.48) 


V6  =  (Ae) 


-N 


h 


V 


f  h  \ 

=  (As)-^ 

-|pc2 

^  -hipu  y 

which  corresponds  to  Xq  =  3.  The  characteristic  solutions  are  normalized  to  satisfy  max  \zk{j)\  =  0(1)  for 

3 

h  tending  to  zero. 

5.4.  Scheme  #  4.  Another  factorizable  scheme  belonging  to  the  family  (4.4)  is  defined  as 

(5.49) 

where  is  similar  to  with  the  desired  discretization  for  the  full-potential  factor  given  as 

(5.50) 

where  is  a  three-point  central  second-order  accurate  approximation  to  the  second  derivative.  The  vector- 
function  qO)  is  defined  similar  to  q^^^ . 

A  particular  solution  can  be  found  in  the  form  of  (5.43)  and  (5.44)  with  the  generalized  symbols 

TH\)=  (u2-c2)a,\(A)  =  %^(i-2  +  A), 

V'^ix)  =  T'^iX)  -  {u‘^d^{X)d<^{X)  -  c^{d%X)f)  . 

The  “maximal-length  footprint”  stencil  of  the  determinant  operator  computed  before  any  cancella- 
tion  occur  includes  seven  points.  Based  on  this  stencil,  the  corresponding  characteristic  equation  is  formed 
as 

(5.52)  u{f3‘-c‘) 


(5.51) 


2h 


0—  -l-  —  —  6^  +  —  10  -l-  3A  -l-  OA^  )  —  0. 


12 


For  the  characteristic  equation,  zero  coefficients  in  the  leftmost  positions  imply  zero  roots,  and,  in  the 
rightmost  positions,  they  imply  infinite  roots.  Six  roots  of  equation  (5.52)  are  Ai,2,3  =  1,  A4  =  |,  A5  =  0, 
and  Ae  =  00. 

The  solution  representation  as  a  linear  combination  of  the  functions  is  relevant  only  for  finite 

Xk  0.  For  Xj.  =  0,  the  corresponding  characteristic  solutions  are  localized  at  the  inflow  boundary,  i.e.,  they 
exhibit  nonzero  values  at  the  inflow  and  are  zero  in  the  interior  and  at  the  outflow  boundary.  By  analogy, 
characteristic  solutions  which  corresponds  to  A^,  =  00  are  localized  at  the  outflow  boundary,  i.e.,  they  are 
nonzero  only  at  some  locations  in  the  vicinity  of  the  outflow  boundary. 

Four  linearly  independent  characteristic  solutions  which  corresponds  to  finite  (nonzero)  Xk  can  be  found 
in  the  usual  form  Zk(j)  =  '^kXk- 


/  jh 

\ 

(5.53) 

Vl  = 

0 

0 

II 

Cl 

> 

,  and  V3  =  p{u^  -  c^) 

1 0  y 

viy 

V  -jhpu 

/ 

correspond  to  Ai,2,3  =  1,  and 


( 

h 

( 

h  \ 

(5.54) 

V4  = 

—hp(?d^{X4) 

1 

= 

[ 

0 

) 

[ 

0  / 

corresponds  to  A4  =  |. 

The  characteristic  solution  Z5  localized  at  the  inflow  (i.e.,  corresponding  to  A5 


0)  is 


(5.55) 


Z5(i) 


f  \ 

\  puh  Sj 


and  the  characteristic  solution  zq  localized  at  the  outflow  (i.e.,  corresponding  to  Xq  =  00)  is 


(5.56) 

where 


Z6(i) 


/  \ 

_  „3li+£  <:N-2 
A'  2  "i 

\  —2>puh  j 


(5.57) 


<5™ 


0,  if  TO  j, 
1,  if  m  =  j. 


Coefficients  Ck,k  =  1, ...  ,6  can  be  found  by  substituting  the  general  solution  into  the  boundary  condition 
equations  that  are  similar  to  (5.42)  with  discretization  "D^  (5.40)  corresponding  to  defined  in  (5.50). 


5.5.  Discretization  Errors.  In  this  section,  the  Loo  norms  of  discretization  errors  in  p  for  Schemes  #  1 
through  #  4  are  compared  for  the  constant-coefficient  problems  corresponding  to  different  Mach  numbers 
(M  =  0.01,0.5,0.99,  and  M*  =  ^  13-^8 ^  0.88).  The  value  of  M*  has  been  chosen  to  illustrate  an  erratic 
convergence  history  for  the  Scheme  #  1.  More  details  are  given  in  Appendix  A.  The  constant  coefficients 
nondimensionalized  with  respect  to  the  density  and  the  speed  of  sound  at  the  sonic  conditions  are  defined  as 


(5.58) 
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Fig.  5.1.  Loo  norms  of  discretization  errors  in  pressure. 


where  7  =  1.4.  The  parameters  of  the  exact  differential  solution  defined  on  the  interval  x  €  [0, 1]  have  been 
chosen  as  =  1,  Cp  =  2,  and  a  =  lir.  The  Loo  norms  of  discretization  errors  in  pressure,  p,  plotted  on 
Figure  5.1  have  been  computed  for  solutions  on  a  sequence  of  uniform  grids  with  h  =  2~‘^,  2“®, . . . ,  2“^®. 

For  low  to  medium  subsonic  Mach  numbers,  the  discretization  errors  of  the  new  factorizable  schemes 
(#  3  and  #  4)  are  very  competitive  with  the  discretization  errors  of  standard  schemes  (#  1  and  #  2)  on 
all  the  grids.  In  this  range  of  Mach  numbers,  the  discretization  errors  of  Scheme  #4  are  smaller  than  the 
discretization  errors  of  Scheme  #3.  The  discretization  accuracy  of  Scheme  #  3  differs  from  the  accuracy  of 
Scheme  #  4  because  of  the  difference  in  the  five-point  and  three-point  approximations  to  the  full-potential 
factor.  The  differences  vanish  as  the  Mach  number  tends  to  unity  because  the  contribution  of  the  full- 
potential  factor  becomes  negligible. 

For  subsonic  Mach  numbers  approaching  unity,  the  characteristic-upwinding  Scheme  #  2  is  obviously 
superior,  exhibiting  discretization  errors  that  are  nearly  two  orders  of  magnitude  smaller  that  the  discretiza¬ 
tion  errors  achieved  by  other  schemes  on  the  same  grids.  Although  not  shown  in  the  figures,  the  situation  is 
similar  for  errors  in  velocity  at  Mach  numbers  approaching  zero.  Scheme  #  2  in  one  dimension  is  very  close 
to  ideal  because  the  entire  system  can  be  cast  as  three  scalar  convection  equations.  A  standard  dimension- 
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by-dimension  extension  to  multiple  dimensions  loses  this  property  as  well  as  the  more  fundamental  property 
of  discrete  factorizability.  Thus,  in  multiple  dimensions,  the  new  factorizable  schemes  proposed  in  this  paper 
are  expected  to  offer  comparable  accuracy  and  considerable  improvements  in  efficiency. 

As  shown  in  Figure  5.1,  convergence  of  discretization  errors  for  Scheme  #  1  may  be  very  erratic.  This 
behavior  is  explained  by  presence  of  the  spurious  global  solutions.  In  particular,  the  choice  of  the  Mach 
number  M  =  M*  corresponds  to  the  spurious  solutions  varying  as  A  pattern  of  discretization  error 

increases  on  each  fourth  grid  is  notable  in  the  test  with  M  =  M* .  Analysis  performed  in  Appendix  A 
confirms  that  the  discretization  error  convergence  rates  for  this  scheme  do  not  settle  to  any  particular  value. 
The  mean  convergence,  however,  obeys  the  order  property.  The  convergence  rate  pattern  is  dictated  by 
the  characteristic  frequency  of  the  spurious  solutions.  The  new  factorizable  schemes  overcome  this  disorder 
and  exhibit  monotonic  convergence  rates  with  an  asymptote  determined  by  the  approximation  order  of  the 
discretizations  in  the  interior. 

Another  interesting  feature  associated  with  the  new  factorizable  schemes  is  the  asymptotic  behavior  of 
the  auxiliary  discrete  function  The  amplitude  of  this  function  is  0{h)  in  some  0(/i)-small  neighborhoods 
of  the  boundaries  and  decreases  exponentially  quickly  to  0(/i^)-size  in  the  interior.  This  behavior  does  not 
compromise  the  second-order  accuracy  in  the  physical  variables  and  It  is  explained  by  interactions 
of  the  interior  discretizations  with  the  overspecified  boundary  conditions.  The  amplitude  of  becomes 
uniformly  0(/i^ )-small  if  another  set  of  boundary  conditions  that  better  suit  the  interior  discretizations  is 
applied. 

6.  Multigrid  Method.  An  FMG  algorithm  solves  the  target-grid  equations  proceeding  from  the  coars¬ 
est  grid  to  finer  grids.  The  goal  of  the  algorithm  is  fast  reduction  of  the  current-grid  algebraic  errors  below 
the  level  of  discretization  errors,  before  interpolating  solutions  to  the  next  finer  grid.  The  algebraic  errors 
on  a  given  grid,  which  are  defined  as  the  differences  between  the  approximate  and  exact  discrete  solutions, 
are  reduced  through  a  Full  Approximation  Scheme  (FAS)  [3,  4,  13,  22,  26,  27]  multigrid  cycle,  in  which 
corrections  to  the  fine-grid  nonlinear  equations  are  obtained  from  coarser  grid  solutions.  The  FMG-FAS 
method  is  described  below  by  means  of  a  two-grid  notation,  in  which  the  fine  grid  is  denoted  by  superscript 
h  and  the  coarse  grid  by  superscript  2h. 

Let  the  fine-grid  nonlinear  equations  be  defined  as 

(6.1)  R‘’(q'‘)  =  f''. 

The  initial  fine-grid  solution  approximation  is  prolonged  from  the  coarse-grid  solutions  as 

(6.2)  q'^=Vq^'^, 

where  V  denotes  a  prolongation  operator  used  by  the  FMG  algorithm  for  solution  interpolation;  the  “hat” 
notation  is  applied  to  distinguish  from  the  prolongation  operator  V  used  within  FAS  multigrid  cycles  for 
interpolation  of  coarse-grid  corrections.  After  forming  the  initial  fine-grid  solution  approximation  q^,  a  two- 
level  FAS  multigrid  cycle  is  applied  as  following:  Several  (or  perhaps  one  or  even  zero)  relaxation  sweeps 
are  applied  to  the  fine-grid  equations  to  obtain  an  improved  solution  approximation  q^.  The  coarse-grid 
equations  at  level  2h  to  be  solved  for  corrections  to  q^  are  defined  as 

(6.3)  R^‘’(q2'‘)  =  R^*’(7^q'‘)  -b  7^(f'‘  -  R‘’(q'‘)), 

where  TZ  and  TZ  denote  restriction  operators  for  transferring  the  fine-grid  solutions  and  residuals  to  the  coarse 
grid,  respectively.  These  coarse-grid  equations  are  then  solved  by  some  iterative  method  (or  directly  if  the 
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grid  is  coarse  enough).  The  corrections  from  the  coarse-grid  (grid  2h)  solutions  are  prolonged  to  the  fine 
grid  as 


(6.4)  q/*  =  q/*  +  p(q2/»  _  JZqhy 

Several  relaxation  sweeps  follow  the  coarse-grid  correction  interpolation  to  complete  one  FAS  multigrid  cycle. 


FV(v„v,)  Cycle 


Fig.  6.1.  Schematic  of  the  f-level  FV(vi,V2)  cycle,  where  vo  denotes  the  number  of  relaxation  sweeps  on  the  coarsest 
mesh  and  denotes  the  grid  with  mesh  spacing  h. 

In  multilevel  versions  of  the  FAS  cycle,  the  coarse-grid  equations  are  themselves  solved  with  7  cycles  of 
the  algorithm  applied  recursively,  where  7=1  would  correspond  to  a  V  cycle  and  7  =  2  to  a  W  cycle.  In  our 
numerical  tests,  a  version  of  the  V  cycle,  termed  FV(i/i,i/2)  cycle,  is  used.  The  multilevel  FV(i/i,i/2)  cycle 
has  been  employed  earlier  in  [24]  and  is  derived  from  the  target-grid  two-level  FAS  cycle  described  above  by 
applying  a  multilevel  FMG  algorithm  with  the  V(i/i,  1/2)  cycle  to  solve  the  coarse-grid  equations.  Parameters 
vi  and  1/2  denote  the  number  of  relaxation  sweeps  on  the  downward  and  upward  legs  of  the  cycle.  The  cycle 
is  sketched  in  Figure  6.1. 

The  distributed  relaxation  method  has  been  applied  as  described  in  Sections  3  and  4.  The  notations  Lh 
and  Mh  are  used  below  for  the  principal  linearization  of  the  operator  R**  of  (6.1)  and  the  corresponding 
distribution  matrix,  respectively.  The  values  of  auxiliary  variables  (5w^  have  been  overspecified  outside  of  the 
computational  domain  by  zeros  wherever  it  is  necessary.  The  convection  operators  in  the  LhMh  matrix  have 
been  solved  by  downstream  marching.  The  full-potential  operator  has  been  relaxed  with  downstream  Gauss- 
Seidel  relaxation  sweeps.  The  inter-grid  transfer  operators  are  the  following:  The  prolongation  operators, 
V  and  V,  are  the  second-order  symmetric  linear  interpolations  of  the  primitive-variable  corrections.  The 
operator  TZ  restricting  residuals  is  the  second-order  full- weighting  operator.  The  solution-restriction  operator 
it  is  the  injection  operator. 

7.  Numerical  Tests.  The  results  of  numerical  solutions  of  nonlinear  nonconservative  equations  corre¬ 
sponding  to  quasi-one-dimensional  subsonic  flow  in  a  convergent-divergent  channel  are  reported  in  this  sec¬ 
tion.  The  differential  equations  are  (5.1)  with  x  €  [0, 1]  and  the  area  distribution  term  a{x)  =  1  —  0.8a;(l  —  a;) 
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(a)  5-pt  full-potential  discretization 

Fig.  7.1.  Maximum  discretization  errors  versus  grid  size. 


Two  factorizable  discrete  schemes  approximating  (5.1)  are  tested 


(7.1) 
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Table  7.1 

The  L2-norm  of  discretization  errors  in  p  and  the  ratio  of  algebraic- to- discretization  errors  for  the  FMG-1  solver  with 
FV(2,2)  cycles. 
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0.01 
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0.2098x10-3 

0.01 
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0.01 

The  principal  linearization  operator  and  the  distribution  matrix  used  in  distributed  relaxation  are  one¬ 
dimensional  versions  of  (4.4)  and  (4.5).  The  discrete  equations  are  solved  with  an  FMG-1  algorithm  em¬ 
ploying  one  FV(2,  2)  cycle  on  each  grid.  The  coarsest  grid  corresponds  to  /i  =  1/16.  The  discrete  boundary 
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conditions  are  overspecified  on  each  grid  from  the  exact  differential  solution  of  the  fully  subsonic  flow  prob¬ 
lem  with  the  inflow  Mach  number  M  =  0.5.  The  maximum  discretization  errors  versus  mesh  size  shown  in 
Figure  7.1  demonstrate  the  second-order  spatial  convergence.  In  Table  7.1,  the  ratios  of  the  algebraic  error 
in  the  L2-norm  to  the  discretization  error  in  this  norm  for  the  FMG-1  algorithm  are  quite  small  at  each  grid. 
The  complexity  of  this  solver  is  about  30  minimal  work  units,  where  a  minimal  work  unit  is  defined  as  the 
number  of  operations  required  for  one  target-grid  residual  evaluation.  The  number  is  relatively  large  as  is 
typical  for  one  dimension;  in  two  dimensions  the  complexity  of  this  algorithm  would  be  about  9.3  minimal 
work  units,  representing  TME  according  to  the  definition  given  at  the  beginning  of  the  paper. 

The  residual  convergence  history  versus  FV(2,2)  multigrid  cycles  are  shown  in  Figure  7.2.  The  equations 
are  numbered  according  to  (7.1).  The  convergence  rates  are  grid-independent  and  are  roughly  one  order  of 
magnitude  per  cycle.  The  discretization  with  the  three-point  full-potential  factor  shows  slightly  better 
convergence  rates,  although,  the  rates  are  slightly  slower  than  the  asymptotic  convergence  rate  expected 
if  solving  only  the  scalar  three-point  full-potential  equation.  Although  not  shown,  the  convergence  rates 
somewhat  deteriorate  when  a  multigrid  V  cycle  is  used  instead  of  the  FV  cycle. 


Fig.  7.2.  Residual  convergence  for  the  FMG  algorithm  with  five  FV(2,2)  multigrid  cycles  at  each  grid  level  for  the 
nonconservative  equations  (7.1). 


8.  Concluding  Remarks.  A  multigrid  method  is  defined  as  having  textbook  multigrid  efficiency 
(TME)  if  solutions  to  the  governing  system  of  equations  are  attained  in  a  computational  work  that  is  a 
small  (less  than  10)  multiple  of  the  operation  count  in  one  target-grid  residual  evaluation.  A  way  to  achieve 
TME  for  the  Navier-Stokes  equations  is  to  apply  the  distributed  relaxation  method  separating  the  elliptic 
and  hyperbolic  partitions  of  the  equations.  Design  of  a  distributed  relaxation  scheme  for  the  Navier-Stokes 
systems  can  be  significantly  simplified  if  the  target  discretization  possesses  two  properties:  (1)  factorizability 
and  (2)  consistent  approximations  for  scalar  factors. 

The  paper  has  introduced  a  new  family  of  factorizable  discrete  schemes  for  the  multidimensional  Eu¬ 
ler  equations.  The  schemes  include  a  mechanism  that  allows  one  to  obtain  any  desired  discretization  for 
the  full-potential  factor  of  the  system  determinant  without  compromising  the  scheme  factorizability.  This 
property  opens  the  door  for  applying  the  distributed  relaxation  technique  leading  to  TME  in  complicated 
computational  fluid  dynamics  simulations. 
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The  accuracy  of  the  new  schemes  has  been  analyzed,  compared  with  the  accuracy  of  standard  schemes, 
and  found  competitive.  TME  and  fast  grid-independent  residual  convergence  rates  have  been  demonstrated 
in  solution  of  fully  subsonic  quasi-one-dimensional  flow  in  a  convergent /divergent  channel. 


Appendix  A.  Erratic  Convergence  of  Basic-Scheme  Discretization  Errors. 

In  this  appendix,  an  analysis  of  discretization  errors  for  the  basic  scheme  (5.20)  is  presented.  This 
analysis  explains  both  the  erratic  convergence  rates  and  the  mean  second-order  convergence  exhibited  by 
the  discretization  errors. 

Taking  the  exact  solution  (5.5)  of  the  differential  problem  as  a  Fourier  component  guarantees  that  the 
particular  solution  (5.24)  is  a  second-order  accurate  approximation,  and  that  the  differences  between  the 
solutions  converge  asymptotically  monotonically.  Thus,  a  deviation  from  the  monotonic  convergence  may 
occur  only  because  of  coefficients  Ck,k  =  1,...,4,  of  the  characteristic  solutions  z^,  to  the  homogeneous 
problem.  For  the  set  of  overspecified  boundary  conditions  (5.29),  the  coefficients  cu  are  found  as 
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and  A3,A4,u(^\  and  are  defined  in  Section  5.1. 

All  the  values  bk  converge  to  zero  with  second-order  rates,  and  this  provides  a  mean  second-order 
convergence  for  discretization  errors.  The  determinant  of  T  is  a  linear  combination  of  some  integer  powers 
(between  —N  and  N)  of  A3.  (Recall,  that  A4  =  AJ^.)  The  sequence  does  not  converge  to  any 

value  as  N  tends  to  infinity  —  rather  it  orbits  the  unit  circle  with  the  frequency  determined  by  4>]  (j),  in 
turn,  is  controlled  by  the  Mach  number.  Therefore,  the  determinant  of  T  is  also  oscillatory.  Theoretically, 
one  could  construct  a  set  of  boundary  conditions  and  choose  the  Mach  number  so  that  T  degenerates  on  a 
particular  grid.  The  discrete  problem  on  this  grid  would  become  ill  posed.  Further  grid  refinement  would 
unavoidably  result  in  repeating  (or  closely  approaching)  ill-posedness  on  an  infinite  number  of  grids.  For 
overspecified  boundary  conditions,  we  did  not  find  a  Mach  number  and  a  grid  to  enforce  the  degeneration  of 
the  matrix  T.  However,  significant  variations  of  the  absolute  value  of  the  determinant  have  been  observed. 

Figure  A.l  presents  a  convergence  history  (circle  symbols)  for  the  Loo  norms  of  discretization  errors  in 
pressure.  The  parameters  of  the  exact  solution  have  been  chosen  as  Cu  =  l,C'p  =  2,  a  =  Ttt.  The  Mach 
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Mach  Number  M  =  M  ~  0.88 


Fig.  A.l.  Loo  norms  of  discretization  errors  (circles)  and  the  history  of  the  determinant  of  T (triangles)  for  the  basic 
discrete  scheme  with  the  overspecified  boundary  conditions. 


number, 


(A.4) 


M*  = 


'  5i2  +  8  +  SVWTl 

25*2  +  16  ’ 


with  t  =  tan(^),  corresponds  to  (f)  =  7r/4.  The  number  of  grid  points  has  been  gradually  increased  (by  1) 
rather  than  doubled  in  passing  to  the  next  fine  grid.  For  comparison,  Figure  A.l  exhibits  the  history  of  the 
determinant  of  T  (triangle  symbols)  as  well.  The  discretization  errors  do  not  converge  monotonically.  Spikes 
of  large  discretization  errors  occur  periodically,  closely  following  the  pattern  of  the  determinant  T  behavior: 
smaller  determinant  values  correspond  to  larger  discretization  errors.  A  mean  second-order  convergence  rate 
is  observed. 
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